rm(list=ls())


# Check that required packages are installed:
want = c("foreign", "ggplot2", "MASS", "reshape2", "plyr", "dplyr", "ggrepel", 
         "lme4", "coda", "relaimpo", "relimp", "arm", "boot", "entropy", "here")
have = want %in% rownames(installed.packages())
if ( any(!have) ) { install.packages( want[!have] ) }
# load packages
junk <- lapply(want, library, character.only = TRUE)
rm(have, want, junk)

options(scipen = 999)

#=============================================================================================================
# Functions
#=============================================================================================================
# Function to extract the variance of random effects from model outputs
extract <- function(x) {
  obj <- VarCorr(x)
  eff <- NULL
  for(n in names(obj)){
    vars <- sapply(obj[n], attr, "stddev")^2
    if(length(vars) > 1){names(vars) <- paste0(n, ".", row.names(vars))}
    eff <- c(eff, vars)
  }
  res <- attr(obj, "sc")^2
  names(res) <- "Residual"
  eff <- c(eff, res)
  data.frame(VarRE = round(eff, 5))
}


#=============================================================================================================
# Load Image
load("allmodels.RData")
#=============================================================================================================

#=============================================================================================================
# Code only correct-not correct
#=============================================================================================================
data$correct_01 <- ifelse(data$correct_p == 1, 1, 0)

#=============================================================================================================
# MODELS
#=============================================================================================================
### Empty model
m.c.logit.empty <- glmer(correct_p_nc ~ 
                           (1|id) + (1|syslab),
                         data = data, family = binomial(link = "logit"))
summary(m.c.logit.empty)
summary(m.c.logit.empty)$logLik
AIC(m.c.logit.empty)
BIC(m.c.logit.empty)

# Store Random effects
ran.eff.m.c <- data.frame(cbind(inv.logit(fixef(m.c.logit.empty)[1] + ranef(m.c.logit.empty)$syslab[, 1]),
                                inv.logit(se.ranef(m.c.logit.empty)$syslab[, 1])))
names(ran.eff.m.c) <- c("COEF", "SE")
ran.eff.m.c$syslab <- rownames(ranef(m.c.logit.empty)$syslab)
ran.eff.m.c <- merge(ran.eff.m.c, ratio.lm.boot)

### Full model
m.c.logit <- glmer(correct_p_nc ~ educ + age_gsd +
                     lrext_cen + polint_cen + polinfo_cen +
                     ch_polar + n_parties_medcen +
                     ri_l2_m +
                     (1|id) + (1|syslab),
                   data = data, family = "binomial",
                   control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1000000)))
summary(m.c.logit)
summary(m.c.logit)$logLik
AIC(m.c.logit)
BIC(m.c.logit)

#=============================================================================================================
# Plot predictions
#=============================================================================================================
n.r <- 100
sim <- data.frame(NULL)
# Save the coefficients from the models
meanbetas <- c(fixef(m.c.logit))
covmatrix <- vcov(m.c.logit)
set.seed(461982)
MCdata <- mvrnorm(n=1000, meanbetas, covmatrix)

md <- 0
partrange <- seq(min(data$ri_l2_m), max(data$ri_l2_m), length.out = n.r)

for (i in 1:n.r){
  pt <- partrange[i]
  x.new <- c(1,1,0,mean(data$lrext_cen,na.rm = T),
             mean(data$polint_cen,na.rm = T),mean(data$polinfo_cen,na.rm = T),
             0,0,pt)
  betahat <- inv.logit(MCdata %*% x.new)
  p <- mean(betahat)  
  lo <- quantile(betahat, probs=c(0.025))
  hi <- quantile(betahat, probs=c(0.975))
  sim <- rbind(sim, c(pt, 2^pt, p, lo, hi))
}
colnames(sim)<-c("X", "expX", "Y", "LO", "HI") 


quartz(type = 'pdf', file = 'FIGURE_E1.pdf',
       width=8, height=7, dpi=300)
ggplot(sim,aes(x = X,y = Y)) + 
  geom_hline(yintercept = mean(sim$Y), color = "gray70", linetype = 1) +
  geom_text_repel(data = ran.eff.m.c, aes(x = ri_l2_m, y = COEF, label = syslab), 
                  size = 3, col = "gray30") +
  geom_line(aes(x=X,y=LO),color="black", linetype = 2) +
  geom_line(aes(x=X,y=HI),color="black", linetype = 2) +
  geom_line(aes(x=X,y=Y),color="black") +
  theme_bw() +
  xlab("Importance ratio (\U03C9)") + 
  ylab("Predicted Probability\n(95% C.I.)")
dev.off()

